home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / roots.i < prev    next >
Text File  |  1996-02-13  |  10KB  |  338 lines
  1. /*
  2.    ROOTS.I
  3.    Collection of root, maximum, and minimum finders.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1996.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. local roots;
  11. /* DOCUMENT roots.i
  12.        defines:
  13.      nraphson     - Newton-Raphson/bisection root solver (scalar)
  14.      f_inverse    - function inverse by Newton-Raphson (vectorized)
  15.      mnbrent      - Brent's method minimizer (scalar)
  16.      mxbrent      - Brent's method maximizer (scalar)
  17.  */
  18.  
  19. /* ------------------------------------------------------------------------ */
  20.  
  21. func nraphson(f_and_dfdx, x0, x1, xerr)
  22. /* DOCUMENT nraphson(f_and_dfdx, x0, x1)
  23.          or nraphson(f_and_dfdx, x0, x1, xerr)
  24.  
  25.      Find a root of a function by Newton-Raphson iteration, backed
  26.      up by bisection if the convergence seems poor.  The subroutine
  27.      F_AND_DFDX must be defined as:
  28.           func F_AND_DFDX (x, &f, &dfdx)
  29.      returning both the function value f(x) and derivative dfdx(x).
  30.      If F_AND_DFDX always returns dfdx==0, nraphson uses bisection.
  31.      The value of x is constrained to lie within the interval from
  32.      X0 to X1; the function values at these two points must have
  33.      opposite sign.  The iteration stops when the root is known to
  34.      within XERR, or to machine precision if XERR is nil or zero.
  35.  
  36.      f_inverse is a "vectorized" version of nraphson.
  37.  
  38.      Based on rtsafe from Press, et. al. Numerical Recipes, Ch 9.
  39.  
  40.    SEE ALSO: mnbrent, mxbrent, f_inverse
  41.  */
  42. {
  43.   if (is_void(xerr) || xerr<0.0) xerr= 0.0;
  44.   x0= double(x0);
  45.   x1= double(x1);
  46.   local f, dfdx;
  47.  
  48.   /* get function value at endpoints -- derivatives unused */
  49.   f_and_dfdx, x0, f, dfdx;
  50.   if (f == 0.0) return x0;
  51.   dxo= f;
  52.   f_and_dfdx, x1, f, dfdx;
  53.   if (f == 0.0) return x1;
  54.   if (f*dxo > 0.0) error, "f(x0) and f(x1) have same sign";
  55.  
  56.   if (f > 0.0) {
  57.     xlo= x0;
  58.     xhi= x1;
  59.   } else {
  60.     xlo= x1;
  61.     xhi= x0;
  62.   }
  63.  
  64.   /* first guess is midpoint */
  65.   x= 0.5*(xhi+xlo);
  66.   dx= dxo= abs(xhi-xlo);
  67.   f_and_dfdx, x, f, dfdx;
  68.   if (f < 0.0)      xlo= x;
  69.   else if (f > 0.0) xhi= x;
  70.   else return x;
  71.  
  72.   for (i=1 ; i<=nr_maxits ; ++i) {
  73.     xo= x;
  74.     if (((x-x1)*dfdx-f)*((x-x0)*dfdx-f) >= 0.0 ||
  75.     abs(2.*f) > abs(dxo*dfdx)) {
  76.       /* take bisection step if N-R step would be out of bounds
  77.      or if previous step did not converge fast enough */
  78.       dxo= dx;
  79.       dx= 0.5*(xhi-xlo);
  80.       x= xlo+dx;
  81.  
  82.     } else {
  83.       /* take N-R step */
  84.       dxo= dx;
  85.       dx= f/double(dfdx);
  86.       x-= dx;
  87.     }
  88.  
  89.     /* quit on either machine precision or requested precision */
  90.     if (x==xo || abs(dx)<xerr) return x;
  91.  
  92.     f_and_dfdx, x, f, dfdx;
  93.     if (f < 0.0) xlo= x;
  94.     else         xhi= x;
  95.   }
  96.  
  97.   error, "nr_maxits iteration count exceeded";
  98. }
  99.  
  100. nr_maxits= 100;
  101.  
  102. func f_inverse(f_and_dfdx, y, x0, x1, xerr)
  103. /* DOCUMENT f_inverse(f_and_dfdx, y, x0, x1, xerr)
  104.          or f_inverse(f_and_dfdx, y, x0, x1, xerr)
  105.  
  106.      Find values of an inverse function by Newton-Raphson iteration,
  107.      backed up by bisection if the convergence seems poor.  The
  108.      subroutine F_AND_DFDX must be defined as:
  109.           func F_AND_DFDX (x, &f, &dfdx)
  110.      returning both the function value f(x) and derivative dfdx(x).
  111.      If the input x is an array, the returned f and dfdx must have
  112.      the same shape as the input x.  If F_AND_DFDX always returns
  113.      zero dfdx, f_inverse will use bisection.
  114.  
  115.      The result x will have the same shape as the input Y values.
  116.  
  117.      The values of x are constrained to lie within the interval from
  118.      X0 to X1; the function value must be on opposite sides of the
  119.      required Y at these interval endpoints.  The iteration stops
  120.      when the root is known to within XERR, or to machine precision
  121.      if XERR is nil or zero.  X0, X1, and XERR may be arrays conformable
  122.      with Y.
  123.  
  124.      f_inverse takes the same number of iterations for every Y value;
  125.      it does not notice that some may have converged before others.
  126.  
  127.    SEE ALSO: nraphson
  128.  */
  129. {
  130.   if (is_void(xerr) || xerr<0.0) xerr= 0.0;
  131.   x0+= 0.0*y;
  132.   x1+= 0.0*y;
  133.   local f, dfdx;
  134.  
  135.   /* get function value at endpoints -- derivatives unused */
  136.   f_and_dfdx, x0, f, dfdx;
  137.   f-= y;
  138.   dxo= f;
  139.   f_and_dfdx, x1, f, dfdx;
  140.   f-= y;
  141.   if (anyof(f*dxo > 0.0)) error, "f(x0)-y and f(x1)-y have same sign";
  142.  
  143.   dfdx= x0;
  144.   mask= double(f > dxo);
  145.   maskc= 1.0-mask;
  146.   xlo= x0*mask + x1*maskc;
  147.   xhi= x1*mask + x0*maskc;
  148.  
  149.   /* first guess is midpoint */
  150.   x= 0.5*(xhi+xlo);
  151.   dx= dxo= abs(xhi-xlo);
  152.   f_and_dfdx, x, f, dfdx;
  153.   f-= y;
  154.   mask= (f < 0.0);
  155.   list= where(mask);
  156.   if (numberof(list)) xlo(list)= x(list);
  157.   list= where(!mask);
  158.   if (numberof(list)) xhi(list)= x(list);
  159.  
  160.   for (i=1 ; i<=nr_maxits ; ++i) {
  161.     xo= x;
  162.     mask= ((((x-x1)*dfdx-f)*((x-x0)*dfdx-f) >= 0.0) |
  163.        (abs(2.*f) > abs(dxo*dfdx)));
  164.     list= where(mask);
  165.     if (numberof(list)) {
  166.       /* take bisection step where N-R step would be out of bounds
  167.      or if previous step did not converge fast enough */
  168.       dxob= dx(list);
  169.       xob= xlo(list);
  170.       dxb= 0.5*(xhi(list)-xob);
  171.       xb= xob+dxb;
  172.     } else {
  173.       xb= dxb= dxob= [];
  174.     }
  175.  
  176.     list= where(!mask);
  177.     if (numberof(list)) {
  178.       /* otherwise take N-R step */
  179.       dxon= dx(list);
  180.       dxn= f(list)/dfdx(list);
  181.       xon= x(list);
  182.       xn= xon-dxn;
  183.     } else {
  184.       xn= dxn= dxon= [];
  185.     }
  186.  
  187.     x= merge(xb, xn, mask);
  188.     dx= merge(dxb, dxn, mask);
  189.     dxo= merge(dxob, dxon, mask);
  190.     /* check for uniform convergence either to requested precision
  191.        or to machine precision */
  192.     if (allof((x==xo) | (dx<xerr))) return x;
  193.  
  194.     f_and_dfdx, x, f, dfdx;
  195.     f-= y;
  196.     mask= (f < 0.0);
  197.     list= where(mask);
  198.     if (numberof(list)) xlo(list)= x(list);
  199.     list= where(!mask);
  200.     if (numberof(list)) xhi(list)= x(list);
  201.   }
  202.  
  203.   error, "nr_maxits iteration count exceeded";
  204. }
  205.  
  206. /* ------------------------------------------------------------------------ */
  207.  
  208. func mxbrent(f, x0, x1, x2, &xmax, xerr)
  209. /* DOCUMENT fmax= mxbrent(f, x0, x1, x2)
  210.          or fmax= mxbrent(f, x0, x1, x2, xmax)
  211.          or fmax= mxbrent(f, x0, x1, x2, xmax, xerr)
  212.  
  213.      returns the maximum of the function F (of a single argument x),
  214.      given three points X0, X1, and X2 such that F(X1) is less than
  215.      either F(X0) or F(X2), and X1 is between X0 and X2.  If the
  216.      XMAX argument is provided, it is set to the x value which
  217.      produced FMAX.  If XERR is supplied, the search stops when
  218.      a fractional error of XERR in x is reached; note that XERR
  219.      smaller than the square root of the machine precision (or
  220.      omitted) will cause convergence to machine precision in FMAX.
  221.  
  222.      The algorithm is Brent's method - a combination of inverse
  223.      parabolic interpolation and golden section search - as adapted
  224.      from Numerical Recipes Ch. 10 (Press, et. al.).
  225.  
  226.    SEE ALSO: mxbrent, nraphson, f_inverse
  227.  */
  228. {
  229.   from_mxbrent= 1;
  230.   return mnbrent(f, x0, x1, x2, xmax, xerr);
  231. }
  232.  
  233. func mnbrent(f, x0, x1, x2, &xmin, xerr)
  234. /* DOCUMENT fmin= mnbrent(f, x0, x1, x2)
  235.          or fmin= mnbrent(f, x0, x1, x2, xmin)
  236.          or fmin= mnbrent(f, x0, x1, x2, xmin, xerr)
  237.  
  238.      returns the minimum of the function F (of a single argument x),
  239.      given three points X0, X1, and X2 such that F(X1) is less than
  240.      either F(X0) or F(X2), and X1 is between X0 and X2.  If the
  241.      XMIN argument is provided, it is set to the x value which
  242.      produced FMIN.  If XERR is supplied, the search stops when
  243.      a fractional error of XERR in x is reached; note that XERR
  244.      smaller than the square root of the machine precision (or
  245.      omitted) will cause convergence to machine precision in FMIN.
  246.  
  247.      The algorithm is Brent's method - a combination of inverse
  248.      parabolic interpolation and golden section search - as adapted
  249.      from Numerical Recipes Ch. 10 (Press, et. al.).
  250.  
  251.    SEE ALSO: mxbrent, nraphson, f_inverse
  252.  */
  253. {
  254.   if (is_void(xerr) || xerr<2.e-8) xerr= 2.e-8;
  255.   while (1.+xerr*xerr == 1.) xerr*= 2.0;
  256.   golden= (sqrt(5.)-1.)/(sqrt(5.)+1.);
  257.   epsilon= 1.e-20;
  258.  
  259.   x0= double(x0);
  260.   x1= double(x1);
  261.   x2= double(x2);
  262.   xlo= min(x0, x2);
  263.   xhi= max(x0, x2);
  264.  
  265.   x= w= v= x1;
  266.   e= 0.0;
  267.   fx= fw= fv= (from_mxbrent? -f(x) : f(x));
  268.  
  269.   for (i=1 ; i<=br_maxits ; ++i) {
  270.     xm= 0.5*(xlo+xhi);
  271.     abserr= xerr*abs(x) + epsilon;
  272.     abserr2= 2.*abserr;
  273.     if (abs(x-xm) <= (abserr2-0.5*(xhi-xlo))) {
  274.       xmin= x;
  275.       return (from_mxbrent? -fx : fx);
  276.     }
  277.  
  278.     if (abs(e) > abserr) {
  279.       /* attempt a trial parabolic fit to (w,v,x) -- the three
  280.          smallest points seen so far */
  281.       r= (x-w)*(fx-fv);
  282.       q= (x-v)*(fx-fw);
  283.       p= (x-v)*q - (x-w)*r;
  284.       q= 2.*(q-r);
  285.       if (q > 0.0) p= -p;
  286.       q= abs(q);
  287.       e1= e;
  288.       e= d;
  289.       if (abs(p)>=abs(0.5*q*e1) || p<=q*(xlo-x) || p>=q*(xhi-x)) {
  290.     /* take golden section step -- parabolic step is crazy */
  291.     e= (x>=xm? xlo-x : xhi-x);
  292.     d= golden*e;
  293.       } else {
  294.     /* take inverse parabolic step */
  295.     d= p/q;
  296.     u= x+d;
  297.     if ((u-xlo)<abserr2 || (xhi-u)<abserr2) d= sign(xm-x)*abserr;
  298.       }
  299.  
  300.     } else {
  301.       /* take golden section step */
  302.       e= (x>=xm? xlo-x : xhi-x);
  303.       d= golden*e;
  304.     }
  305.  
  306.     /* evaluate at next point, avoiding meaninglessly small step size */
  307.     u= x + (abs(d)>=abserr? d : sign(d)*abserr);
  308.     fu= (from_mxbrent? -f(u) : f(u));
  309.  
  310.     if (fu <= fx) {
  311.       /* new point is best so far */
  312.       if (u >= x) xlo= x;
  313.       else        xhi= x;
  314.       v= w;   fv= fw;
  315.       w= x;   fw= fx;
  316.       x= u;   fx= fu;
  317.  
  318.     } else {
  319.       if (u < x) xlo= u;
  320.       else       xhi= u;
  321.       if (fu<=fw || w==x) {
  322.     /* new point is second best so far */
  323.     v= w;   fv= fw;
  324.     w= u;   fw= fu;
  325.       } else if (fu<=fw || v==x || v==w) {
  326.     /* new point is third best so far */
  327.     v= u;   fv= fu;
  328.       }
  329.     }
  330.   }
  331.  
  332.   error, "br_maxits iteration count exceeded";
  333. }
  334.  
  335. br_maxits= 100;
  336.  
  337. /* ------------------------------------------------------------------------ */
  338.